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Abstract. In this paper I give a detailed account of an ab initio methodology for 
describing strong electronic correlations in nanoscale devices hosting transition metal 
atoms with open d- or /-shells. The method combines Kohn-Sham Density Functional 
Theory for treating the weakly interacting electrons on a static mean-field level with 
non-perturbative many-body methods for the strongly interacting electrons in the open 
d- and /-shells. An effective description of the strongly interacting electrons in terms 
of a multi-orbital Anderson impurity model is obtained by projection onto the strongly 
correlated subspace properly taking into account the non-orthogonality of the atomic 
basis set. A special focus lies on the ab initio calculation of the effective screened 
interaction matrix U for the Anderson model. Solution of the effective Anderson model 
with the One-Crossing approximation or other impurity solver techniques yields the 
dynamic correlations within the strongly correlated subspace giving rise e.g. to the 
Kondo effect. As an example the method is applied to the case of a Co adatom on the 
Cu(OOl) surface. The calculated low-bias tunnel spectra show Fano-Kondo lineshapes 
similar to those measured in experiments. The exact shape of the Fano-Kondo feature 
as well as its width depend quite strongly on the filling of the Co Sd-shell. Although this 
somewhat hampers accurate quantitative predictions regarding lineshapes and Kondo 
temperatures, the overall physical situation can be predicted quite reliably. 


1. Introduction 

Modern experimental techniqnes now allow to reliably create, manipnlate and control 
nanoscale devices with atomic precision in the lab thns bringing the dream of molecular 
electronics or nanoelectronics to create ultimately miniaturized electronic devices from 
single molecules closer to reality m [21 m m. Prospective building blocks for molecular 
electronic circuits such as molecular rectihers m and held-effect transistors have 
already been demonstrated in experiments. The use of magnetic atoms or molecules 
promises to further enhance the functionality of molecular devices by exploiting the 
spin-degree of freedom of the electron in addition to its charge. Such devices could 
serve e.g. as basic building blocks for nanoscale spintronics applications [a [IQ] or as 
ultimately miniaturized magnetic information storage devices [H]. 

Naturally, quantum effects play a crucial role in electronic devices of such tiny 
dimensions. Consequently, experiments with atomic- and molecular-scale devices have 
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produced a wealth of quantum phenomena such as conductance quantization [T2] . 
quantum interference [I31II3, or quantum phase transitions |T5]. On the other hand, 
details of the atomic structure also play an important role for determining the electronic 
properties of nanoscale devices, especially regarding the contact between molecule and 
metal leads [ini [m HE]. Also the coupling to the leads can significantly alter the 
electronic and magnetic properties of nanoscale devices by broadening and shifting of 
energy levels, as well as screening effects. Hence a proper theoretical description of 
nanoelectronic devices needs to take into account all of the following: quantum effects, 
the actual atomic structure of the device and the coupling to the leads. 

The now standard approach for the description of molecular electronic devices is 
to combine density functional theory (DFT) calculations with the Landauer transport 
theory or with the non-equilibrium Green’s function formalism (NEGF) [T^ 1^ 1^1^ . 
The DFT based transport approach yields an effective mean-field description for the 
electronic structure and transport properties of molecular devices, taking into account 
quantum effects, as well as the actual atomic structure of the device, and the coupling 
of the device to the metallic leads. The approach works quite well for the description of 
metallic nanocontacts and nanowires and carbon nanotubes [211 ESI EE]. On the other 
hand, it was realized quite early on that this approach often overestimates conductances 
of molecules attached to metal leads by orders of magnitude. Its origin has been a matter 
of debate for over a decade and is still not completely settled [211 EEl EE] • 

Moreover, nanoscale devices comprising magnetic atoms or molecules often display 
phenomena induced by so-called strong dynamic correlations that arise when the 
effective Goulomb interaction between the electrons exceeds their kinetic energies. 
Dynamic correlations can have a profound impact on the electronic and magnetic 
structure and the transport properties of the system. One of the most intriguing 
phenomena induced by dynamic correlations in nanoscale devices is probably the Kondo 
effect [23 EH]: Below a critical temperature characteristic of the system, the Kondo 
temperature T^, the atomic or molecular spin forms a many-body singlet state with 
the nearby conduction electrons, thereby screening the magnetic moment of the device. 
The correlations usually originate from the strongly interacting open Sd- or 4/-shells 
of transition metal atoms. But also molecular orbitals of purely organic molecules only 
weakly coupled to the leads can give rise to strong correlations. This is corroborated 
by the fact that the Kondo effect is not only frequently observed in molecular devices 
comprising transition metal atoms [29] [301 ISIl [321 ESI EH EEl EEl E3 ESI EHl HD] , but also 
for devices made from purely organic molecules [HIEll 1121 SSIIHIIIE]- By construction 
the DFT based transport method being a static mean-field approach cannot c^ture 
the dynamic correlations that lead e.g. to the Kondo effect in nanoscale devices [j This 

f Recently it has been shown by Bergfield et al. that the exact exchange correlation functional yields 
the exact transmission at the Fermi level in the case of the simple Anderson impurity model. However, 
even the exact Kohn-Sham spectrum does not yield a correct description of the spectral function and 
transmission outside the Fermi level. Hence the renormalization of the Kondo peak by the interactions 
cannot be captured by Kohn-Sham DFT based transport calculations |46) . 
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neglect of dynamic correlations could also be behind the afore mentioned overestimate 
of the conductances of molecular devices by the DFT based transport approach since 
dynamic correlations can lead to a strong renormalization of the quasi particles relevant 
for the transport through the molecule [17]. 

Recent efforts to go beyond the DFT based transport approach are to combine 
time-dependent DFT (TDDFT) with the NEGF [HI 09] or the GW approximation 
with NEGF [501 ED E2]- A problem of the TDDFT approach is that the standard 
approximations for TDDFT functionals in connection with an adiabatic exchange 
correlation kernel does not yield an improvement for the description of correlation 
effects with respect to the static DFT approach. Some progress has been made 
recently in that direction by Ending a non-adiabatic exchange-correlation kernel for 
strongly correlated systems but only in the context of simplified models such as the 
Hubbard or Anderson model [HSllM]. The GW based transport approach on the other 
hand has been implemented in a fully ab initio way and has been applied to realistic 
molecular devices. Although GW yields an energy-dependent self-energy for describing 
the electronic interactions and thus captures dynamic correlation effects to some extend, 
it is perturbative in nature and thus strong electronic correlations such as those leading 
to the Kondo effect or the Mott-Hubbard metal-insulator transition are not properly 
described. 

Here I give a detailed account of a different ab initio approach for the description of 
strongly correlated molecular conductors which has been developed, successively refined 
and extended in previous work [55l |56l ED ED EH]- fhis approach only the strongly 
interacting part of the electronic spectrum is described by advanced many-body methods 
in order to capture dynamic correlations effects. The weakly to moderately interacting 
part of the electronic system is still treated on a static mean-field level by standard Kohn- 
Sham DFT (KSDFT). This approach is basically an adaption of the DFT-|-Dynamical 
Mean-Field Theory (DFT-I-DMFT) approach [HD ED ED E3]) which has been developed 
for the realistic description of strongly correlated solids, to the special situation of 
nanoscale conductors. Similar approaches for treating strong correlations in molecular 
devices have recently appeared in the literature [HD ED ED ED ED EH] . 

This paper is organized as follows: In Sec. [D a detailed account of the so far 
developed methodology is presented. In Sec. [3] the methodology is applied to the 
case of a Go adatom at the Gu(OOl) surface which has been studied extensively in 
the recent past [ZD [ZD CD ES] and thus presents an ideal testbed for the theory. In 
Sec. 01 I draw conclusions from comparison of the results to the experiments and other 
theoretical methods. I also discuss some of the caveats of the developed theory and 
possible solutions to these problems as well as future directions. 

2. Methodology 

The typical situations encountered in experiments with atomic and molecular devices 
are depicted in Fig. [D (a) A magnetic molecule suspended between the tips of a metal 
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(b) 



Figure 1. Typical situations encountered in molecular electronics/spintronics: (a) 
A magnetic molecule bridging the tips of a nanocontact, (b) A magnetic molecule on 
a metal surface probed by an STM tip. (c) Schematic sketch of model that captures 
both situations shown in (a) and (b). A central atom or molecule (turquoise) hosting 
strongly correlated levels C (red) is connected to two metal leads L and R (yellow). 
The device region D (blue) contains the central atom/molecule and part of the leads. 
The polarization region P (magenta) extends over that part of the atom/molecule and 
the lead(s) in close proximity to C. 


nanocontact and (b) a magnetic atom or molecule deposited on a metallic surface 
probed by an STM tip. The magnetism and hence the strong correlations of the 
molecule are here assumed to stem from a single transition metal atom at its center. 
But the approach can be easily generalized to the case of multiple magnetic atoms 
by adaption of the Dynamical Mean-Field Theory (DMFT) to the case of molecular 
conductors |S7|. Both situations depicted in Figs. [I](a,b) can be described by the model 
depicted schematically in Fig. [11(c): the central region, called device region D, contains 
the molecule or atom and part of the leads, and will be described on the level of KSDFT. 
Within the atom or molecule the correlated subspace C yields the strongly interacting 
levels of the atom/molecule that will be treated by advanced many-body techniques 
in order to capture the strong dynamic correlations. The polarization region P where 
the polarizability is calculated in order to compute the screened interaction U of the 
strongly correlated subspace on the other hand extends over that part of the molecule 
and/or leads in immediate vicinity of the correlated subspace C. This approach can in 
principle also be applied directly to the case of purely organic molecules. In that case 
one has to identify the molecular orbitals responsible for the strong correlations [74] . 

The approach has been implemented within the ANT.G package m which 
interfaces the Gaussian quantum chemistry code [76] in order to implement the DFT 
based ab initio transport methodology for molecular conductors. The Gaussian code 
makes use of Gaussian atomic orbitals as basis sets for performing quantum chemistry 
and DFT calculations of finite clusters and molecules. The ANT.G package embeds the 
hnite cluster representing the device region into bulk electrodes in order to model the 
transport situation depicted in Fig. [H However, the formalism developed below is not 
specihc to Gaussian basis sets. It can directly be applied to any atomic basis set, as 
for example the Fireball orbitals used in the SIESTA code mi. Even more general, the 
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formalism might be applied to any basis set as long as the different snbspaces (D,P and 
C) can be dehned in a meaningfnl way. 

2.1. Non-orthogonal basis sets and projection onto a subspace 

We now have to carefnlly dehne the projections onto the different snbspaces taking 
into acconnt the non-orthogonality of the atomic basis set. The choice of projection 
strongly inflnences physical qnantities associated with the snbspace snch as the density 
and electronic occnpancy of the snbspace as has been shown recently by Soriano and 
Palacios pS]. 

We assnme that the Hilbert space H of onr system is spanned by a (hnite) set of 
non-orthogonal orbitals H = {|a)}, i.e. H = span(if), and {a \ (3) = Sag 7 ^ 0 for 
I a), 1/9) G H. We now want to project onto a snbspace M of H spanned by a snbset 
M = {|m)} of the orbitals |a) G H, i.e. M G H. Dne to the non-orthogonality of the 
orbitals |a) G iP, snbspace M will in general have a hnite overlap with the snbspace R 
spanned by the rest of the orbitals |r) G i? = H\M, i.e. Smr = I 7^ 0 1^) ^ ^ 

and |r) G R. Hence the qnestion arises how to dehne a proper projection Pm onto that 
snbspace. We note that there has actnally been some controversy abont this qnestion 
in the literature (see e.g. Ref. ua and references therein). 

It turns out that the proper choice for Pm is actually quite obvious: Let us hrst 
consider the simplest case of the subspace M being spanned by a single orbital \m). 
By dehnition, the projection operator for a single state is simply P^ = |m)(m|. This 
dehnition is independent of how (in which basis) the Hilbert space of the entire system 
is dehned; i.e. it does not matter whether \m) forms part of the basis set spanning 
the entire Hilbert space or not; or in case it does whether it has some overlap with the 
Hilbert space R spanned by the rest of the basis set. 

Hence it is clear that the projection -Pm for the subspace M can be written 
in an orthonormal basis set M-*- = {I’ti-*-)} spanning the subspace M as 

-Pm = Such an orthonormal set can always be found by 

Lowdin orthogonalization of the original non-orthogonal set spanning M: \m^) = 
where Sm is the overlap matrix between the basis set elements of 
M only and is an abbreviation for (Sm)”^'^^, i-e. the matrix power —1/2 of the 

matrix Sm- Hence we hnd for the projection operator: 

Pm = Y. Y. H (Sm ^^)„imx(SM ^^)^x„|m)(n| 

m,n&M 

= l"i)(SM )mn(ri| (1) 

m^n£M 

which is nothing but the identity operator for the subspace M written in the non- 
orthogonal basis set. It has been argued on more formal grounds that this choice for the 
projection is actually the only physical reasonable one as it is the only one that leads 
to a tensorial consistent occupancy matrix which generates a Hermitian potential pg. 
Note that the subspace projection Pm dehned here corresponds to the projector with 
regard to the A metric denoted by Pm [ZH]- 
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Also note that in general we cannot write the identity operator for the entire system 
as the sum of the projection onto subspace M and subspace R spanned by the rest of 
the basis set if there is some overlap between the two subspaces, i.e. I ^ Pm + Pr- 
Rather we have to correct for the overlap between the two subspaces: 

i= Y. |«)(S-'W(/5| = Pm + Pr + 0 (2) 

a,/ 3 £H 

where is the inverse of the overlap matrix for the entire system. O is an operator 
correcting the sum of projections by the overlap between the two subspaces M and 
R. = Pr + O dehnes the projection onto a new subspace M which is actually 
orthogonal to subspace M. The projection Pm thus dehnes an orthogonalization scheme 
which orthogonalizes R with respect to subspace M preserving the latter. 

Now let us have a look at how an operator A acting on the full Hilbert space is 
projected onto the subspace M: 

-4m = PmAPm = Y \m) {m'\A\n') {Sfl)n'n{n\ 

= Y |^^)(Sm AmSm )mn(n| = Y \m){AM)mn{n\ (3) 

m,nEM m,n€M 


where Am = {{m\A\n)) is the direct matrix given by the matrix elements of A with the 
basis {|q;)} of subspace M, and Am = SJ^^AmS^^ is the so-called nuclear matrix in that 
basis. Note that for an orthonormal basis of M we have Am = Am- 

Frequently, we will also have to project an operator A given for some subspace M 
onto a smaller subspace M' C M; 


Am'= Pm'AmPm'= Y PM'\m) (Am) ayinlPw 

m,n€M 

= Y I m){AM)mn{n I q'){S^,)q>n'{n'\ 


m,nGM 
m’ 


— Y I"^ 0 (Sm'Sm'm Am SmM' 


(4) 


where Sm'm is the overlap matrix between orbitals \m') G M' and orbitals |m) G M 
Hence we obtain the following expression for the nuclear matrix of subspace M' in terms 
of the nuclear matrix for subspace M: 

Am' = Sjyj/ Sm'M Am SmM' Sjyj/ (5) 


On the other hand, we may also have the opposite situation where we have some 
operator Am only dehned on subspace M, and we want to know the direct matrix for 
the entire space H, i.e. 

{a\Au\P) = XI I m){Au)mn{n \ ft) (6) 

m,nEM 

Hence the direct matrix of the operator Am is given by 


Am — ShmAmSmh 


(7) 
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The central quantities both in DFT based transport calculations of molecular electronics 
devices and in quantum many-body theory are Green’s functions (GF). The one-body 
GF is dehned as the resolvent of the one-body Schrodinger equation |80j : 

G{z){z + fi-H)=i (8) 

where z is complex, /i is the chemical potential, and H is the Hamiltonian of the system. 
G{z) has poles at the eigen values of H for a hnite system or a branch cut on the real 
axis at the energy bands for an inhnite system. Its spectral representation in terms of 
the eigen states \k) oi H is given by: 

G(z) = (z + ^,-Hr^ = J:Af2A_ (9) 

^ Z + fX — Ck 

The GF operator projected onto subspace M is given by: 

= PmG{z)Pm = |Q:)(GM(2:))a/3(/d| (10) 

a,l3£M 

Dehning the GF of the isolated subspace M as 

guiz) = {{z + fi)pM — Hm) ^ (11) 

and the self-energy operator Em associated with the coupling of the subspace M to the 
rest of the world as 

Sm(^) = [9m{z)] ^ — [Gm(^)] ^ (12) 

it is possible to rewrite the projected GF as 

Gm{z) = (^{z + P)Pm — Hm — Em(^)) (13) 

The self-energy Em(^) is not to be confused with the one describing electron-electron 
interactions in the many-body GF formalism. Note that in many-body physics in the 
context of the Anderson impurity model [ 8 T] Sm(^) is often called hybridization function 
and is denoted by Am(: 2 ). 

One can easily write Em(^) in terms of the GF for the isolated (i.e. not coupled 
to M) complementary space M dehned by Pjq, Qm^z) = {{z + fi)Pu — -^m)”^ ns 
Sm(^) = -^m,m^m( 2 :)-^M,M where Hm,m = PmHPm = (Pm,m)^- In order to hnd the 
matrix representations Gm and Gm of the projected GF Gm( 2 :), eq. (B, is multiplied 
with the denominator of the r.h.s., the matrix elements are taken and the subspace 
identity Pm is inserted between the two factors of the l.h.s.: 

{alGuiz) I«')(Sm )a/3'(/9'l {{z + /i)PM - Sm(;^)) \f3) = {a\pM\^) 

(14) 

Hence we hnd for the direct GF matrix 

Gm(^) = Sm {{z + /i)SM — Hm — Sm(^)) ^ Sm 

^ {{z + - Hm - Sm(^)) 


(15) 
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and for the corresponding nnclear matrix 

Gm(^) = (( 2 : +/i)SM — Hm — Sm(^)) ^ (16) 


2.3. Many-body Green’s functions and Feynman diagrams in an atomic basis set 


The generalization of the one-body Green’s fnnction for an (effectively) non-interacting 
system to the case of interacting electrons are the single-particle Green’s fnnction or 
single-particle propagators. The single-particle Matsnbara GF [H2] for atomic states a 
and a' is dehned as 




-{Tr 


[Cc 


r 


(r')]) 


(17) 


where r is imaginary time and the creation and annihilation operators obey the 
generalized anti-commntation rnles for non-orthogonal basis sets [83] : 


■(cq., Oq )■ S^ 


af) 


(18) 


The Fonrier transform with respect to imaginary time r yields the Matsnbara GF for 
imaginary freqnencies (called Matsnbara freqnencies): 

rd 

Gaafi^) = / dr Q,/(r, 0) (19) 

Jo 

By analytic continnation to the real freqnency axis one obtains the retarded single¬ 
particle GF uj + ir]). 

Gaa> dehnes the direct single-particle GF matrix G. In the absence of interactions 
the single-particle GF matrix G tnrns ont to be eqnal to the one-body GF matrix 
dehned as the resolvent of the one-body Schrodinger eqnation. Analogons to the one- 
body GF we can also dehne the nnclear matrix for the interacting single-particle GF as 
G = S-^GS-F 


For the development of a diagrammatic expansion for the interacting GF in a 
non-orthogonal basis in terms of the Gonlomb interaction and the non-interacting GF, 
one has to either use the nuclear GF matrix in combination with the direct Gonlomb 
interaction matrix, or the direct GF matrix in combination with the nuclear matrix of 
the Gonlomb interaction [83]. Here we will work with the nuclear matrix for the Green’s 
functions and the direct matrix for the interactions. 

The bare Gonlomb interaction in an atomic basis set is given by: 


ye-e 


^a.P\ct'l3' C, 


t J 

OLG^Ol' g' 




( 20 ) 


where Vag-^a'y is the nuclear matrix of the Gonlomb interaction [H3], i.e. V(l,2) = 
S(1)“^S(2)“^V(1, 2)S(2)“^S(1)“^ and the direct matrix elements are given by 


14 


a(3;a'^‘ 


f = e 


||ri - r2\\ 


( 21 ) 


The Feynman diagrams for the GF and the Gonlomb interaction in an atomic basis 
set are shown in Fig. |2l 
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Figure 2. Feynman diagrams for the single-particle Green’s function G and the bare 
Coulomb interaction 14 in an atomic basis set. 


2.4- DFT based transport calculations 

We consider the situation schematically depicted in Fig. [Tt^c). The central device region 
D containing a molecule is coupled to two electrodes L and R. This situation can be 
realized in a number of ways as shown in Figs. [^a,b): (a) A molecule bridging the 
tips of a nanocontact or (b) a molecule deposited on a metal substrate and coupled to 
an STM tip. In addition to the molecule the device region D contains those parts of 
the two electrodes which are in close proximity to the molecule and whose electronic 
structure is modified by the presence of the molecule and vice versa. In the case of the 
molecular bridge (a) the tips of the nanocontact are included in the device region while 
in the case of the molecule on the substrate (b), part of the surface and of the STM tip 
are included in the device region. 

The electronic structure of the central device region is calculated ab initio on the 
level of DFT in the Kohn-Sham (KS) framework, taking into account the coupling to 
the electrodes L and R. The Kohn-Sham Green’s function of the device region D is given 
by: 

= {{z + /i)SD - H° - Sl(z) - Sr(z))-' (22) 

where is the KS Hamiltonian of the device region which yields an effective mean-held 
description of the electronic structure of the device region. Sl(^) and Sr( 2 :) are the lead 
self-energies associated with the coupling of the device region to the bulk electrodes. 

From the device GF the electronic density can easily be calculated by integration 
up to cj = 0 (corresponding to the chemical potential /i) 

Dq = —Im— [ dujG^{u + ir]) (23) 

TT J—oo 

the density matrix yields a new KS Hamiltonian for the device region thus closing the 
self-consistency cycle of the KS calculation. Hence we can self-consistently calculate 
the electronic structure of the device region taking into account the coupling to the 
electrodes (open system). 

In contrast to D, the electronic structure (Hamiltonian) of the electrodes L and 
R, and hence the self-energies are kept hxed during the self-consistent calculation of 
the electronic structure of D. Depending on the situation, different models for the 
bulk electrodes can be employed. One can for example choose nanowires |H1], embed 
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the cluster into a perfect crystalline surface calculated ab initio [^, or use so-called 
absorbing boundary conditions (ABC) [85]. Here we choose a tight-binding Bethe lattice 
model [Ml with realistic tight-binding parameters obtained from DFT calculations [HZ] • 
The actual choice of the electrode model is not crucial for calculations as long as the 
bulk electrodes are far enough away from the central scattering region, i.e. the device 
region is chosen big enough and contains a sufficiently big part of the electrodes [88]. 

Once the KS calculation is converged the transport properties can be calculated 
within the Landauer approach from the transmission function which is given by: 


T\u) = Tr[rL(a;)G°^(a;)rR(a;)G° (a;)] 


(24) 


where Fl = ^(Sl — Sj)) and Fr = ^(Sr — Sr) are the so-called coupling matrices which 
yield the broadening of the device region due to the coupling to the leads. From the 
transmission function the current and conductance can be calculated using the Landauer 
formula 

I { y ) = ^ y dujT°{uj) (/(a; - /xr) - /(a; - /xr)) (25) 

where /il and /xr are the electrochemical potentials of the left and right lead, respectively, 
dehned by the applied bias voltage eV = /xl — hn- Note that in general the transmission 
function T^{oj) also depends on the applied voltage V, i.e. = T^{u,V), and 
actually has to be calculated out of equilibrium by combining the KSDFT with the 
NEGF [ini EH] El] • However, within the mean-held like KSDFT based NEGF approach 
the transmission is often not so strongly voltage dependent, and hence current and 
conductance can be approximated well by the equilibrium transmission T°(a;, 0) at least 
for sufficiently small bias voltages. 

In the typical situation of an STM setup (Fig.[T](b)), most of the applied bias voltage 
V will drop near the sharp STM tip, i.e. the electrochemical potential of the substrate 
remains hxed to the equilibrium one figuh = h while that of the STM tip changes with 
the bias /xtip = + eV. The differential conductance for low bias at zero temperature is 

then directly given by the transmission function: 


G(I/) = 




2e d 


reV Op2 

/ dujT^{uj) = — xT^{eV) 
Jo h 


(26) 


In contrast, for the situation of a molecule coupled symmetrically to two leads (Fig.[T](a)), 
the voltage will drop more or less symmetrically across the junction, i.e. = fi — eV 12 
and /iR = yU -|- eV 12. Hence for the conductance we obtain now 


G(V) = I X 


d /•+eG2 


dujT^\ 


.2 r 


iV) = 






+ T° - 


eK 


(27) 


dV J-eV/2 h 

In a more general situation where the coupling is neither completely symmetric nor 
completely asymmetric, more sophisticated modelling of the electrostatics or even a 
KS-NEGF calculation would be necessary in order to hnd the actual voltage drop. 
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Next we have to identify the strongly correlated subspace C. Usually C will be formed 
by the open d- or /-shells of transition metal atoms. However, also molecular orbitals of 
purely organic molecules such as Ceo or carbon nanotubes weakly coupled to electrodes 
can show strong correlations if the effective interaction in these levels is big in comparison 
with the broadening due to the coupling to the leads. Our approach is completely general 
in this respect. 

From now on we assume that the orbitals 0 forming the subspace C are mutually 
orthogonal (but not necessarily orthogonal to the other orbitals in the device region). 
This can always be achieved by simple Lowdin orthogonalization of subspace C. 
However, note that often the orbitals spanning C are already mutually orthogonal. For 
example in the case of the atomic orbitals forming the open d- or /-shell of a transition 
metal atom, or in the case of molecular orbitals which are the eigenstates of the KS 
Hamiltonian of the molecule and hence by construction are orthogonal. In order to 
account for the strong correlations in subspace C an effective Coulomb interaction term 



e 


^ ijki 

n-n^ 


(28) 


is added acting on the orbitals in C. Note that Uik-ji is not the bare Coulomb interaction 
but an effective interaction which is usually much lower than the bare one due to 
screening processes by the conduction electrons. In the next section it is shown how 
to calculate Uik-ji ab initio from the DFT electronic structure. The full many-body 
Hamiltonian of the strongly interacting subspace C now reads: 


7/c = (29) 

where the one-body part Hq = J2ij,a{4>i\HQ\(j)j)cl„Cj^ is given by projection of the KS 
Hamiltonian onto C. However, since the Coulomb interaction has been taken into 
account already on a mean-held level in the Kohn-Sham Hamiltonian, we also need to 
subtract a double-counting correction (DCC) term: 

= PcH^Pc - (30) 

Unfortunately, the DCC term is not exactly known for DFT, and several 
approximation schemes are used in practice |H9]. Here the so-called atomic limit or 
fully localized limit (FLL) is employed [90] , but generalized to the case of an anisotropic 
Coulomb repulsion Uu-jj [59] : 

(v^h = E Uuu, ("i - ^) - (31) 

where Uj = (c]cj) is the electronic occupation of orbital 0^, Mq is the dimension of 
subspace C, Jh is the Hund’s rule coupling given by the orbital-averaged exchange 
matrix elements Uij.ji, and Nq = J2jec '^j is the total electronic occupation of subspace 

C. 
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According to fflB]) the self-energy (a.k.a. the hybridization fnnction) associated with 
the conpling of C to the rest of the system is given by 

Ac(a;) = (a; + /i)lc - - [G° (a;)]"' (32) 

where the projected GF of the correlated snbspace C C D can be calculated from the 
device GF according to ([5]) as 

Gc(i^) = ScdGd(w)Sdc (33) 

As customary in many-body physics we will call Ac(cj) the hybridization function from 
now on. The many-body Hamiltonian 'He of subspace G together with the hybridization 
function Ac(ci;) dehne a multi-orbital Anderson impurity model (AIM). Solution of the 
AIM yields the self-energy Sc(ci;) describing the strong electronic correlations within 
the G subspace which is fed-back to the DFT calculation in order to obtain electronic 
spectra and transport properties of the molecular device (see Sec. 12.8|) . 

3.6. Computation of the effective interaction in the correlated subspace 

The effective interaction between the electrons in the correlated subspace G is not 
the bare Goulomb interaction because of screening processes by formation of electron- 
hole (e-h) pairs in the rest of the system. Therefore the screened Goulomb matrix 
elements Uik-ji are considerably lower than the bare Goulomb interaction Vik-ji- The 
screening of the bare interaction by formation of e-h pairs can be calculated within the 
so-called Random Phase Approximation (RPA) (see e.g. the book by Mahan [H2] or any 
other textbook on quantum many-body theory). However, screening of the electrons 
within the G subspace will already be taken into account by the impurity solver. Hence 
the contribution of the impurity subspace G to the screening needs to be subtracted 
out. By doing so one arrives at the so-called constrained Random Phase Approximation 
(cRPA) |9I]. 

In order to calculate the effective screened interaction Uik-ji of subspace G within 
cRPA we hrst dehne the so-called polarizability region P in which screening processes 
due to formation of e-h pairs are taken into account for calculating the screened 
interaction. P comprises the strongly correlated subspace G and a sufficient portion 
of the surrounding atoms of subspace G as is schematically indicated in Fig. [H In 
principle, the whole D region could be chosen as P. However, in practice this is often 
not feasible because of computational limitations if the device region is reasonably big. 
Also as it turns out the screening of the correlated subspace G by the surrounding 
conduction electrons is relatively localized due to the usually localized nature of the 
strongly correlated orbitals making up G. 

Within RPA the screened interaction W is given by the Dyson equation shown 
diagrammatically in Fig. |3] which in an atomic orbital basis set can be written 
algebraically as 

hFQ:i/3i;02/32 (H ) Tz) Va-ijSy,a2fS2 ^ *^(D D) 
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Figure 3. Dyson equation for RPA screened interaction for atomic basis set. Wiggly 
lines correspond to the bare Coulomb interaction V, double wiggly lines to the RPA 
screened interaction W. 
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^lUi^2^2 


Vr 




dr (Hp) 




(^1: ^ E 2 ) 


(34) 


For the screening of the bare Conlomb interaction V only screening processes within 
region P are taken into acconnt. Hence we have to calcnlate the polarizability (i.e. the 
bnbble diagram in Fig. E]) projected onto the P region: 




(35) 


where the projected GF for the P region can be obtained from the device GF according 
to (jSj) as 


Gp = S;‘SpDGDSDpSf‘ 


(36) 


For a stationary Hamiltonian we can replace the two times in the screened 
interaction and polarizability by time differences: n(ri, T 2 ) —>■ n(ri—r 2 ) and kF(ri, T 2 ) —)■ 
W (ti — T 2 ). Hence (by setting r 2 = 0 and after some renaming), we can write the Dyson 
eqnation for the RPA screened interaction as: 


hFai/3i;a2d2 ^i/3i;«2d2 ^ ^('^) 

rP 


1111^11121^2 


dr' (Hp) 


llltli-,112112 




)^tl2ll2',0!2l32 ) 


(37) 


Here we will only consider the static limit of the screened interaction, i.e. 
pj/o ^ W{uj = 0) = f dr W{t). Becanse of the /3-periodicity of n(r) we also have 
Jq drn(r — r') = Jq drU^r) = H^. Hence we obtain the following Dyson eqnation for 
the static screened interaction hF°: 


aiPi;a2P2 


^Oiil3i-,a202 A E n 

ll\Ul^2l'2 


ai0i-,^iiviijP'P ) r-Oil\V'2'i2 ^l2U2\a2 02 


0^ lyO 


(38) 


The static Polarizability H^ is now fonnd easily by integrating a Green’s fnnction 
prodnct over the freqnency domain: 


(np)^ll/i;^2l'2 — / dr {Yi-p, 


tliuv,(12112 


v)= / <iTy:(G?);„,(-T)(G“p 


0 \(J 




j) 


' iu 


ILOn ^ 

1 

27r J—00 


duj E(G°p);».(*‘^) (Gp);.„(“) 


(39) 
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Dyson equation for fully screened RPA interaction Wq of subspace C in 


terms of effective interaction U. Orbital indexes have been suppressed here. 


where in the last step we have taken the zero temperature limit [jd — )■ cxd) rendering the 
discrete Matsubara frequencies continuous. 

We now dehne superindices / := (* 1 ,^ 2 ) in order to rewrite the Dyson equation in 
form of a matrix equation. Hence we have W° = {Wj j) etc., and the Dyson equation 
can be written in matrix form as: 

w° = v + vn°w° (40) 

Solving for the static screened interaction W° we hnd: 

W°=(l-Vn°)"V (41) 

Projection to the correlated subspace C then yields the RPA screened interaction 
for the correlated electrons Wq. However, since the screening within the correlated 
subspace will already be taken into account by the impurity solver in a more or less 
exact way, the screening of the correlated electrons by themselves has to be subtracted 
out in order to obtain the effective interaction U. Hence the effective interaction U 
is the partially screened interaction that results in the fully RPA screened interaction 
Wq when taking into account only the polarizability He within the C subspace. The 
corresponding Dyson equation is shown diagrammatically in Fig. 01 Solving for the 
effective interaction U we arrive at the following “unscreening” equation [92l |93l [M] 
computing U from W^: 

U = W° (lc + n°W°)-i (42) 


In order to determine the screening within subspace C, we have to calculate the 
polarizability corresponding to subspace C. 

__ 1 roo ____ __ 

n"c = ^/ (43) 

ZTT J—00 ^ 

where Gc(*cu) = ScdGd(^c<;)Sdc- _ _ 

It is important to realize that 11^ and Gc are not just submatrices of the 
corresponding bigger matrices lip and Gp in the P region of the device due to the 
overlap between the subspaces. Neglecting this detail can result in serious errors in the 
computation of the effective Coulomb interaction U; Due to the numerical instability 
of eq. (142]) small inaccuracies in computing 11^ can result in large errors and even in 
completely unphysical effective interactions. The numerical instability of eq. (021) can 
be seen by rewriting it as 


U 




n -1 


+ n« 


-1 


(44) 
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As the fully screened interaction Wq is usually quite small (compared to the bare 
Coulomb interaction) and positive, [W^]“^ is big and positive. On the other hand the 
screening of the correlated electrons by themselves is usually strong, and therefore IIq is 
big and negative. Hence in order to obtain U we are basically subtracting two relatively 
big numbers and inverting the resulting small number so that relatively small errors in 
calculating W° or nc can result in quite large errors for U. It should be noted here 
that in the case of a semiconducting or insulating substrate or host material, as well 
as in the case of insulating compounds this issue is less problematic since then at low 
energies around the Fermi level, the two subspaces are completely decoupled, leading 
to weaker “self-screening” of the correlated electrons, and hence smaller numbers for 
[W^]-^ and n^. 

However, in the case of a metallic host or substrate considered here, it is thus 
crucial to correctly perform the projections of the different quantities involved in the 
calculation of the effective interaction (Green’s functions, polarizability) to the P and 
C subspaces in order to reliably calculate U. Also the usual neglect of certain product 
basis states in the computation of the screened interaction |95] might be problematic 
in this context. One way to stabilize the numerical evaluation of (H2ll is to decouple 
the correlated subspace from the rest of the system both in the calculation of W and 
of U as proposed by Miyake et al. [96]. This way the self-screening of the correlated 
electrons is reduced considerably, leading to smaller values of and 11°, and 

thus enhancing the numerical stability. However, this can lead to far too high matrix 
elements for the direct Coulomb interaction as will be shown in Sec. El Apparently, 
“mixed propagators” between the correlated subspace and the rest of the system (which 
vanish when the subspaces are decoupled) can be quite important for the screening of 
the effective interaction. 

2.1. Solution of the Anderson impurity model: One-Crossing Approximation 

Since the interaction Uijki is strong in comparison with the single-particle broadening 
(given by the imaginary part of Ac(c(;)), the AIM problem cannot be solved by standard 
perturbation theory in the Coulomb interaction. Instead more advanced many-body 
methods usually starting from an exact diagonalization of the full impurity Hamiltonian 
Hq have to be employed in order to properly take into account the strong correlations 
within subspace C. Here I use the One-Crossing Approximation (OCA) [97] which is 
an improvement over the Non-Crossing Approximation (NCA) [98119^ flOO] . However, 
it should be emphasized that the methodology presented so far can in principle be 
combined with any other method for solving the AIM, as e.g. continuous time Quantum 
Monte-Carlo (CTQMC) [lOlj . or numerical renormalization group (NRG) [102] . or the 
Lanczos diagonalization scheme [66]. 

One advantage of OCA over other schemes is that spectral data can be calculated 
directly on the real frequency axis. Hence in contrast to the numerically exact 
CTQMC, for example, it does not suffer from artifacts introduced by numerical analytic 
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continuation of the spectra from the Matsubara axis to the real axis. Also spurious 
features in the spectra coming from the approximation of the inhnite and continuous 
conduction electron bath in the Anderson model by a hnite and discrete one as in direct 
diagonalization schemes such as Lanczos, are not a problem for OCA since the bath 
is not truncated or discretized. On the other hand, in contrast to the basically exact 
but computationally very demanding NRG, OCA can actually be applied to realistic 
Anderson models of Sd- and 4/-impurities with 5 and 7 impurity-levels, respectively. 

However, being an approximate method, OCA also suffers from some dehciencies 
that one should be aware of. First, as in the case of the simpler NCA, spurious non- 
Fermi liquid behaviour is obtained in the zero-temperature limit, resulting in artifacts 
in the spectral density for low temperatures. While in NCA these artifacts already 
appear below T/^, in OCA the critical temperature below which the artifacts appear is 
signihcantly lower (1-2 orders below Tk)- Another problem of NCA and OCA is the 
violation of certain sum rules especially in the case of multi-orbital Anderson models 
that lead to errors in the high frequency expansion of the electronic self-energy [103] . 
Again, these errors are much less pronounced in OCA than in NCA. 

The basic idea of both NCA and OCA methods is to treat the coupling of the 
correlated subspace C to the rest of the system given by the hybridization function 
Ac(c(;) as a perturbation to the dynamics within the subspace induced by the strong 
electron-electron interactions which is treated exactly. Hence the starting point is an 
exact diagonalization of the many-body Hamiltonian of the correlated subspace: 

Uc = ^Em\m){m\ (45) 

m 

where \m) are the many-body eigenstates of "He and Em the corresponding eigen- 
energies. 

It is now convenient to represent the many-body eigenstates of |m), in terms of 
auxiliary helds or pseudo-particles (PPs) dm, hm which obey (anti-)commutation rules 
depending on the number of electrons represented by the corresponding many-body 
state \m). The physical electron operators Qo-, cI^. are related to the PP operators by: 


= ( 46 ) 

m,n 

where R™”' are the matrix elements of the electron annihilation operator with the many- 
body eigenstates of C: F™” = {m\cia\n). Since the PPs obey (anti-)commutation rules 
a diagrammatic expansion of PP propagators in terms of the coupling to the rest of the 
system is possible. The full PP propagator corresponding to a many-body state |m) is 
then given by 


where Flm(cu) is the PP self-energy describing the dynamic interaction with the other 
PPs induced by the hybridization with the rest of the system (bath). 

NCA consists in an inhnite resummation of self-energy diagrams where conduction 
electron lines do not cross (hence the name). These are the diagrams shown in the left 
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Figure 5. Diagrams for pseudo-particle self-energies in NCA and OCA 
approximations, full lines correspond to conduction electron propagators coupled to 
impurity levels a, double dashed lines to full pseudo-particle propagators. 


box in Fig. [5] for a certain PP m representing a many-body state of N electrons. The 
NCA diagrams describe processes where a single electron (hole) jumps from the bath 
to subspace C and back thereby temporarily creating a PP with N-l-1 (N-1) electrons. 
Hence the NCA self-energy is given by a convolution of the hybridization function Ac(ci;) 
with the PP propagators Gm' of the PPs m' coupled to m. OCA additionally takes into 
account diagrams where two bath electron lines cross as shown in the right box of Fig. |5l 
The algebraic expressions for the OCA self-energies involve double convolutions of two 
hybridization functions with three PP propagators. The exact algebraic expressions can 
be found in the literature [971 [93]. Since the self-energy of a PP m depends on the 
dressed propagators of the other PPs m' that interact via Vhyb with m, the NCA/OCA 
equations have to be solved self-consistently. 

Once the NCA/OCA equations have been solved, the real electronic quantities can 
be calculated from the PP propagators by expanding the real electron operators in terms 
of PP operators by (Hbll . Within NCA, the real electron spectral function is obtained 
from the PP spectral functions as 

pM = 7^ S / An-(0J + e) (48) 

\^ / mm' 

where Am{uj) = is the PP spectral function for PP m and Q is the PP 

charge which is obtained by integration of the PP spectral functions. Again, in OCA the 
expression for calculating the electronic density Pia{oo) is more complicated, involving 
double convolutions of PP spectral functions. From the electron spectral density Pia{^) 
being the imaginary part (modulo tt) of the electron Green’s function Gjo-(i^) we can 
calculate the real part of Giaioj) by Kramers-Kronig. Finally, from the GF Gci.^) 
the electronic self-energy describing the dynamic correlations within G is obtained by 
Sc(w) = [G^(a;)]“^ — [Gc(a;)]“^ where ~ + h‘)Pc ~ Hq — h.Q{u))~^ is the bare 

propagator of subspace G. For a more detailed account of the NGA, OGA and other 
methods based on a hybridization expansion of atomic states see e.g. Refs. [6T], [93] . 
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2.8. Feedback of the self-energy: correlated electronic structure and transport properties 


Once we have solved the Anderson impurity model for the strongly interacting subspace 
C coupled to the rest of the system, we obtain the electronic self-energy describing the 
strong dynamic correlations within subspace C: 


= E B Pc(i*<)l., 01 


(49) 


!jec 


Note that Sc = Sc since we have assumed the basis to be orthonormal within subspace 
C. This self-energy is now fed back to the DFT part in order to obtain the correlated 
electronic structure and transport properties of the system. More specihcally, we obtain 
the correlated device GF 

-1 




(60) 


where is the DCC operator which like Sc('^) only acts on C. According to ([7]) the 
corresponding nuclear matrix of the device GF is 

, -1 


GD(a;) = 


V^IS 


CD 


(51) 


([G»nH|-'-SDclSc(i*<)- -c. 
where the overlap matrices Sdc and Sen sandwiching Sc(a;) — account for the 
overlap between the correlated subspace G and the rest of the system (see eq. ED- 

From the correlated device GF Gd(ci;) we can calculate the correlated electronic 
density analogously to fl2^ by integration of Gd(co’) up to 0 energy: 


Dd = -Im- / 

TT J —oo 


du GD(to’ -|- ip) 


(52) 


From the correlated density in turn a new KS Hamiltonian for the device region can 
be calculated, from which a new correlated density is obtained and so forth until self- 
consistency is reached. Hence we can calculate the effect of the correlation within the 
G subspace onto the charge distribution of the device region. This part corresponds to 
the so-called “charge self-consistency” loop within the DFT-I-DMFT scheme [63] . 

Following Meir-Wingreen [104] . the low-bias transport properties can be obtained 
in complete analogy to the case of KS-DFT transport eqs. (1241127p even in the presence 
of strong correlations from the correlated transmission function 

T{uj) = Tr[rL(a;)GUa;)rR(a;)GD(a;)] (53) 


Note that the strong correlations giving rise e.g. to the Hondo effect are actually 
contained in T{u) via the correlated GF GD(tj). 

In the next section we will see that the Fano-Kondo lineshapes measured by STM 
spectroscopy of magnetic atoms and molecules on metal substrates can indeed be 
reproduced by calculating the conductance from the (zero-bias) transmission function. 
This is due to the fact that the Hondo effect is a low-energy phenomenon, i.e. the Hondo 
peak is observed for very small bias voltages so that hnite-bias effects only play a minor 
role. For the description of actual non-equilibrium phenomena the formalism has to be 
generalized to include the effect of hnite bias voltages. As shown by Meir and Wingreen 
in their landmark papers [104] this can be achieved by generalization of the formalism 
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Nd 

^2 

xz 

yz 

x2 — J/2 

xy 

Sd 

LSDA 

8.13 

1.59 

1.66 

1.66 

1.33 

1.89 

0.82 

LDA 

8.25 

1.66 

1.63 

1.63 

1.52 

1.81 

- 

OGA 

8.26 

1.34 

1.94 

1.94 

1.08 

1.97 

0.86 

OGA (+0.4eV) 

8.15 

1.11 

1.96 

1.96 

1.06 

1.97 

0.94 


Table 1. Total and orbital resolved occupations and spin of Co Sd-shell within DFT 
on the level of LSDA and LDA and DFT+OCA calculations. In the last line we show 
the DFT+OCA results for the Co Sd-levels ed shifted by 0.4eV upwards in energy with 
respect to the FLL. 


to the Keldysh contour. However, in this case the Anderson impurity problem has to 
be solved out of equilibrium which is computationally extremely demanding. So far it 
has only been achieved in the context of the single-level AIM [105[ I106[ 1107] . but not 
for realistic cases. 

3. Results: Co adatom at the Cu(OOl) surface 

Now the developed methodology is applied to the case of a Co adatom deposited on 
the Cu(OOl) surface. This system is an ideal testbed for the theory as it has been 
measured extensively in the recent past [7D1 1108[ 17^ I109[ 1731 IHOj - Fig. [6]^a) shows the 
atomic structure of the device region. The device contains the Co atom on three layers 
of the Cu(OOl) surface and an STM tip consisting of a Cu pyramid grown in the (001) 
direction. The Co atom and its four nearest neighbour Cu atoms have been relaxed with 
GaussianOO m using the local spin density approximation (LSDA) and the LANL2DZ 
double-zeta valence plus outer core electron basis set with core pseudo potentials |lllj 
while the rest of the device atoms have been kept hxed. The interlayer and intralayer 
distances for the hxed Cu atoms are those of a perfect Cu surface taken from Ref. |112j . 
In good agreement with Ref. |113[ 1110] , I hnd that the Co atom relaxes at a height of 
about l.sA above the four nearest neighbour Cu atoms while these in turn are pushed 
by O.IA into the substrate. 

Using ANT.G and the LANL2MB minimal basis set including valence and outer 
core electrons with pseudo potentials m the electronic and magnetic structure 
structure of the device coupled to the tip and substrate electrodes is calculated within 
DFT on the level of LSDA. The Go atom is essentially in a conhguration with 

the two holes in the 3d-shell giving rise to an a approximate spin-1 state of the Go 
atom (see Tab. [T]) again in good agreement with |113] . LSDA basically predicts a mixed 
valence situation for all the orbitals with the individual occupations around 1.6 with 
the exception of the x|/-orbital which is nearly full. 

From the LSDA electronic structure the effective Goulomb interaction Uij-ki for 
the Go 3(i-shell is calculated as described in Sec. EJl For the P region we take into 
account substrate atoms up to the 3rd nearest neighbour, i.e. the 9 Gu atoms closest 
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density-density interaction (eV) 
z"^ xz yz x^ — y"^ xy 

Hund’s coupling (eV) 
z"^ xz yz — y"^ 


5.38 4.27 4.27 3.45 3.46 


xz 

4.27 5.56 3.86 3.73 3.74 

0.60 

yz 

4.27 3.86 5.56 3.73 3.74 

0.60 0.83 

— 2/^ 

3.45 3.73 3.73 5.23 4.28 

0.94 0.82 0.82 

xy 

3.46 3.74 3.74 4.28 5.26 

0.92 0.83 0.83 0.48 


Table 2. Direct Coulomb repulsion matrix elements Uu-^kk (density-density 
interaction) and exchange matrix elements Uik-^ki (Hund’s rule coupling) of effective 
Coulomb interaction for Co 3d-shell. 


to the Co adatom. The change in U from taking into acconnt 2nd nearest neighbonrs 
to 3rd nearest neighbours is about 2%. For 4th nearest neighbours (18 atoms in total) 
the super matrices in the RPA equation flTT|) become too big (linear matrix dimension 
234^ = 54756) to be handled. 

Tab. [2] shows the matrix elements of the effective Coulomb interaction, namely 
the direct Coulomb repulsion matrix elements (density-density interaction) Uu-kk and 
the exchange interaction matrix elements (Hund’s rule coupling) Uik-ki- The average 
density-density interaction is f7 = 4.14 eV. It is strongly screened by the conduction 
electrons, resulting in a reduction of over 80% compared to the bare value of 22.9 eV for 
the Co 3(i-shell. On the other hand, the Hund’s rule coupling is much less affected by 
the screening: it is only reduced by about 10% from its bare value of 0.85 eV, resulting 
in an average Hund’s coupling of Jh = 0.77 eV. Note that the inter-orbital Coulomb 
repulsion {Ua-kk for i ^ k) is related to the average intra-orbital Coulomb repulsion for 
both orbitals and the Hund’s rule coupling via Uu-kk = {Uu-u -b Ukk-,kk)/‘^ — ‘^Uik-ki- 

Both density-density interactions and Hund’s rule coupling are somewhat 
anisotropic (i.e. orbital-dependent). The intra-orbital Coulomb repulsion Uu-u deviates 
only by up to 0.17eV or by to 3% from its mean value of (7 = 5.4eV. The variation is 
stronger for the inter-orbital Coulomb repulsion Uu-^kk, deviating by up to 0.43eV or by 
up to 11% from its mean value of U' = 3.85eV. The Hund’s rule coupling Uik-,ki deviates 
even stronger by up to 0.29eV or by up to 38% from its mean value of Jh = 0.77eH. It is 
worth noting at this point that the complete decoupling of the correlated subspace from 
the rest of the system as proposed in Ref. [96] in order to achieve a stable computation 
of the effective interaction in the case of “entangled bands” produces a much higher 
density-density interaction of about 12eV. Apparently the screening effects of “mixed 
propagators” between the correlated subspace and the rest of the system are actually 
quite important and cannot be neglected. 

Next, the electronic structure of the system is calculated for the paramagnetic case 
on the level of the local density-approximation (LDA) in order to obtain the KS energy 
levels of the Co 3(i-shell e° and hybridization functions Ad{u}) in the absence of spin- 
polarization. Fig. [6])b) shows the imaginary parts of the hybridization functions Arf(a;) 
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Figure 6. Results for Co adatom at Cu(OOl) surface, (a) Atomic structure of 
device part. The Co adatom is shown in grey, (b) Orbitally resolved imaginary part of 
hybridization functions for Co Sd-shell. (c) Orbitally resolved OCA spectral functions 
for Co 3d-shell at T ^ lOK. (d) Total occupation of Co 3d-shell as a function of energy 
shift Ae. (e) Half-width of Kondo feature in spectral function as a function of 
the total shift Ae of Co 3d-levels with respect to energy levels given by FLL DCC. 
(f) Spectral functions of Co z^-orbital for different energy shifts Ae at T ^ lOK. (g) 
DFT-I-OCA transmission functions for different energy shifts Ae (line colours as in 

(f))- 
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for each of the Co 3(i-orbitals which yields the (dynamic) broadening of the orbitals dne 
to the conpling to the snbstrate. We see that the broadening near the Fermi level is 
basically featnreless indicating conpling to the delocalized Cn 4s-states of the snbstrate. 
As can be seen the degenerate xz- and |/z-levels conple most strongly to the these states. 
Becanse of their shape these two orbitals conple very well to the 4s-states of the fonr Cn 
atoms directly nnderneath the Co adatom. On the other hand, the conpling of the xy- 
orbital to the snbstrate is the weakest for all hve orbitals since the direct conpling to the 
nnderneath Cn atoms is strongly snppressed dne to symmetry reasons. The conpling of 
the and the — |/^-orbitals to the snbstrate is intermediate between these two cases. 
At negative energies, the conpling to the localized Cn Sd-states of the snbstrate leads 
to strong peaks in the hybridization fnnctions at energies between -5 and -2 eV. Less 
prononnced peaks at positive energies above 4eV indicate conpling to the Cn 4p-orbitals 
of the snbstrate. 

The bare energies of the Co 3d-levels constitnting the impnrity shell in the 
Anderson impnrity model are obtained from their KS energies = PdH^Pd corrected 
by a DCC term, as explained earlier in Sec. 12.51 The so-called FLL generalized to an 
anisotropic (i.e. orbital-dependent) density-density interaction Uu-^kk is employed (l3lll . 
The valnes for the direct Conlomb repnlsion are the ones shown in Tab.[2l For the Hnnd’s 
rnle conpling the orbital averaged exchange interaction is taken, i.e. Jh = 0.77eV. 

The Anderson impnrity problem presented by the interacting Co 3(i-shell conpled 
to the snbstrate is now solved within OCA as described in Sec. 12.71 For the effective 
Conlomb interaction of the Co 3d-shell we take into acconnt the density-density 
interactions Uu-kk as well as the exchange interactions Uik-ki as given in Tab. [2l At the 
energy levels for the Co 3d-orbitals given by the FLL-DCC fl3Tl) . the total occnpancy 
for the Co 3d-shell is abont 8.2 electrons similar to the ones of the LDA and LSDA 
calcnlations (see Tab.[T]). However, the individnal occnpancies of the 3(i-orbitals are now 
qnite different from the DFT ones, namely they are now closer to integer occnpancies, as 
opposed to the mixed-valence sitnations obtained in the DFT calcnlations. In particnlar, 
the — y^-orbital is now basically half-hlled, and the xz-, yz- and xy-orbitals are nearly 
fnll now. The 2 ;^-orbital is now also closer to half-hlling than before bnt still has strong 
charge flnctnations (occnpancy~ 1.3). Similar to LSDA, the spin of the Co 3(i-shell is 
fonnd to be Sd ~0.87, close to a spin-1 conhgnration. 

In Fig. [6](c) the caicniated spectral fnnctions of the Co 3(i-orbitaIs Pd{u;) (at 
T ~ 1077) are shown. We can see a very strong Kondo peak at the Fermi level in 
the ^^-orbital. The npper Hnbbard peak is here qnite close to the Kondo peak at the 
Fermi level dne to the strong charge flnctnations. This orbital is still qnite close to a 
mixed-valence sitnation. The x^ — |/^-orbital despite being half-hlled and thus bearing a 
spin-1/2 does not yield a Kondo peak. We are dealing here essentially with a so-called 
underscreened Kondo effect |11411115( 1116] : Despite the relatively similar hybridization 
of the z^-channel and the x^ — j/^-channel, the Kondo temperature of the z‘^- 

channel is much higher than that of the x^ — |/^-channel, due to its stronger 

charge huctuations. Hence at hnite temperature T with Tx^x^-y^ < T < only the 
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spin-1/2 in the z^-channel is Kondo-screened, while the spin-1/2 in the channel 

remains unscreened. 

The half-width of the Kondo peak is about 90 K, in very good agreement with the 
experimentally observed values [701 [721 [73]. However, the lineshape of the calculated 
transmission function (red curve in Fig. [6](g)) is rather peak-like, different from the 
experimentally observed asymmetric Fano-lineshaphes. As the DCC for DFT is not 
exactly known and eq. f[3T|) is only an approximation, we now shift the Co 3d-levels 
upwards in energy by an amount Ae thus emptying the Co 3(i-shell as can be seen 
in Fig. [6](d). Emptying the Co 3(i-shell mainly lowers the occupancy of the z^-orbital 
reducing the charge fluctuations for that orbital, while the occupancies of the other 
orbitals are quite stable. Fig. [6](f) shows the effect of shifting the Co 3(i-levels and 
the concomitant reduction of charge fluctuations on the the spectral function of the 
z^-orbital (at T ~ lOK): As the z^-orbital is emptied, its occupation approaches 1, the 
Kondo peak becomes smaller, and the upper Hubbard peak moves away from the Fermi 
level. The width of the Kondo peak decreases at hrst and then starts to grow again for 
shifts > 0.2eV, as can be seen in Fig. [6](e). Note that the non-monotonic behaviour of 
the width of the Kondo peak is actually a hnite temperature effect: As the actual Kondo 
temperature decreases with decreasing charge fluctuations, the Kondo peak in the hnite 
temperature spectra (here T ~ lOK) does not attain its full (zero-temperature) height 
anymore. As the height of the (hnite-T) Kondo peak decreases, its width starts to grow 
again at some point. Hence the half-width of the Kondo peak measured at some hnite 
temperature really only yields an apparent Kondo temperature. 

Fig.[6](g) shows the ehect of shifting the Co 3(i-levels on the low-energy transmission 
spectra. As said above, for the 3(i-levels at the values given by the FLL-DCC the 
transmission function near zero energy is rather peak-like, unlike the ones observed 
experimentally. But when shifting the 3d-levels upwards in energy the lineshapes become 
more asymmetric Fano-like. Good agreement between theoretical and experimental 
Fano-lineshapes is achieved for shifts between 0.4eV to O.SeV. In this regime the 
half-width of the Kondo peak is between 67K and 86K, in good agreement with 
the experimentally observed ones between 70K and lOOK for the Co on Cu(OOl) 
system [TO] [TS] [73] . These results are quite different from those obtained recently for the 
case of Co on Cu(lll) with a similar approach |117] where all the Co 3(i-orbitals give 
rise to Kondo-like resonances at the Fermi level. The reason could be the altogether 
quite different geometric situation at the (111) surface leading to a decidedly different 
symmetry and occupancy for the Co 3(i-shell. 

4. Conclusions 

In conclusion, an ab initio methodology has been developed for describing the impact 
of strong electronic correlations on the electronic structure and transport properties of 
nanoscale devices. Starting from the DFT electronic structure of an embedded nanoscale 
device, an Anderson impurity model is constructed by projection of the Kohn-Sham 
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Hamiltonian onto the correlated subspace. The effective Coulomb interaction U for the 
correlated subspace (impurity) is calculated ab initio from the DFT electronic structure 
by making use of the constrained RPA approach. The solution of the Anderson impurity 
model yields the dynamic correlations originating from strong interactions within the 
correlated subspace in form of a self-energy which is fed back to the DFT calculation in 
order to obtain the correlated electronic structure and transport properties. 

The methodology has been tested for the case of a single Co adatom on Cu(OOl) 
substrate. On a qualitative level the results are in good agreement with experiments: 
A Fano-Kondo feature with the width in good agreement with experiments is obtained 
in the calculated low-energy tunnelling spectra. However, the lineshape of the Fano- 
Kondo feature is not correctly reproduced at the energies for the Co Sd-levels given by 
the double-counting correction. Only when shifting the Co Sd-levels slightly upwards 
in energy good agreement with the experimentally observed lineshapes is achieved. It 
is a well known problem of DFT-I-U and DFT-I-DMFT approaches that the double¬ 
counting correction is not exactly known and in general does not yield the correct 
position (and thus charge) of the correlated levels. Nevertheless, the so-called fully- 
localized limit employed here, is actually not too far off as only moderate shifts are 
necessary to achieve good quantitative agreement with experiments. Importantly, the 
physics is actually not affected by the shifting of the Co Sd-levels: Independent of the 
shift (in that energy range) the Co 3d-shell constitutes essentially a spin-1 system that 
experiences an underscreened Kondo effect. The shifting only affects the weight of the 
Kondo peak by lowering the charge fluctuations in the Kondo-screened orbital. 

Hence the developed methodology is capable of qualitative predictions of strong 
correlation phenomena. But accurate quantitative predictions for example of Kondo 
temperatures and the exact shapes of Fano-Kondo features are difficult as these are 
dependent on the exact occupancy of the correlated subspace which cannot be calculated 
accurately because of the approximate nature of the double-counting correction in our 
approach. One possibility to overcome these difficulties is to make use of the GW 
approach instead of DFT for the description of the weakly interacting part of the system, 
similar to the GW-I-DMFT approach for strongly correlated materials |118[ I119( 1120] 
since for GW the double-counting correction term is exactly known. 
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